!********************
MODULE metric
!********************
USE ElliDef
USE NodeInfoDef
USE libinterp
IMPLICIT NONE
SAVE
PRIVATE
REAl (KIND=qPrec) :: deltax,deltay,deltaz,dmax,xmin,max,ymin,ymax,zmin,zmax
INTEGER nx,ny,nz

REAL (KIND=qPrec), DIMENSION (:,:), POINTER :: h,hx,hy,hxx,hyy

PUBLIC g,g_mask,ind2cord,setprob,DEFINE_TOPO



CONTAINS

SUBROUTINE setprob(Info)

TYPE (NodeInfo), TARGET :: Info

! set global variables in this module from values in Info

deltax=Info%dX(1);deltay=Info%dX(2);deltaz=Info%DX(3)
dmax=Info%dmax
nx=Info%mX(1);ny=Info%mX(2);nz=Info%mX(3)
xmin=Info%xLower(1);ymin=Info%xLower(2);zmin=Info%xLower(3)
IF(.NOT.ASSOCIATED(Info%topo)) ALLOCATE(Info%topo)
! set the local pointers to the topography in Info 
h=>Info%topo%h;hx=>Info%topo%hx;hy=>Info%topo%hy;
hxx=>Info%topo%hxx;hyy=>Info%topo%hyy;



END SUBROUTINE setprob

SUBROUTINE DEFINE_TOPO(Info,time)

TYPE (NodeInfo) :: Info
REAL (KIND=qPrec), OPTIONAL, INTENT (OUT) :: time
REAL (KIND=qPrec) :: dummy
REAl (KIND=qPRec), DIMENSION(:), ALLOCATABLE :: x,y,xg,yg
REAL (KIND=qPrec), DIMENSION (:), ALLOCATABLE :: w
REAL (KIND=qPrec), DIMENSION (:,:), ALLOCATABLE :: depth

INTEGER, DIMENSION (:), ALLOCATABLE :: iw
INTEGER :: lw, liw, intpol(2)=(/3,3/),ier,nxg,nyg,i,j

IF(PRESENT(time)) THEN
CALL CPU_TIME(dummy)
CALL CPU_TIME(dummy)
END IF

! define local coordinates
ALLOCATE(X(0:Info%mX(1)+1))
ALLOCATE(Y(0:Info%mX(2)+1))

x=Info%Xlower(1)+((/0:Info%mX(1)+1/)-.5d0)*deltax
y=Info%Xlower(2)+((/0:Info%mX(2)+1/)-.5d0)*deltay


! load topography
OPEN(99,file='topo.data',status='old',form='unformatted')
READ(99) dmax;Info%dmax=dmax
PRINT*,'max depth =',dmax
READ(99) nxg,nyg

ALLOCATE(xg(nxg),yg(nyg),depth(nxg,nyg))
READ(99) xg
READ(99) yg
READ(99) ((depth(i,j),i=1,nxg),j=1,nyg)


CLOSE(99)

! check coordinates with topography
xg=xg/dmax;yg=yg/dmax;depth=depth/dmax

IF(minval(x)<minval(xg)) PRINT*,'x < xg'
IF(maxval(x)>maxval(xg)) PRINT*,'x > xg',maxval(x),maxval(xg)
IF(minval(y)<minval(yg)) PRINT*,'y < yg'
IF(maxval(y)>maxval(yg)) PRINT*,'y > yg'

! interpolate topography
IF(.NOT.ASSOCIATED(Info%topo)) ALLOCATE(Info%topo)
ALLOCATE(Info%topo%h(0:Info%mX(1)+1,0:Info%mX(2)+1),Info%topo%hx(0:Info%mX(1)+1,0:Info%mX(2)+1),Info%topo%hxx(0:Info%mX(1)+1,0:Info%mX(2)+1))
ALLOCATE(Info%topo%hy(0:Info%mX(1)+1,0:Info%mX(2)+1),Info%topo%hyy(0:Info%mX(1)+1,0:Info%mX(2)+1))


h=>Info%topo%h;hx=>Info%topo%hx;hy=>Info%topo%hy;
hxx=>Info%topo%hxx;hyy=>Info%topo%hyy;


lw=4*(2*SIZE(X)+SIZE(Y))
ALLOCATE(w(lw))
liw=SIZE(X)+SIZE(Y)
ALLOCATE(iw(liw))

CALL RGRD2(nxg,nyg,xg,yg,depth,SIZE(X),SIZE(Y),X,Y,h,intpol,w,lw,iw,liw,ier)


IF(ier/=0) THEN 
   PRINT*,'Error ',ier,' in Interpolating topography'
   STOP
ENDIF

! calculate derivatives
hx=zero
hx(1:Info%mX(1),0:Info%mX(2)+1)=(h(2:Info%mX(1)+1,0:Info%mX(2)+1)-h(0:Info%mX(1)-2,0:Info%mX(2)+1))/(2*Info%dX(1))


hx(0,:)=(h(1,:)-h(0,:))/(Info%dX(1));
hx(Info%mX(1)+1,:)=(h(Info%mX(1)+1,:)-h(Info%mX(1),:))/(Info%dX(1))
hy=zero
hy(0:Info%mX(1)+1,1:Info%mX(2))=(h(0:Info%mX(1)+1,2:Info%mX(2)+1)-h(0:Info%mX(1)+1,0:Info%mX(2)-1))/(2*Info%dX(2))


hy(:,0)=(h(:,1)-h(:,0))/Info%dX(2)
hy(:,Info%mX(2)+1)=(h(:,Info%mX(2)+1)-h(:,Info%mX(2)))/Info%dX(2)





DEALLOCATE(x,y,xg,yg,w,iw,depth)

IF(PRESENT(time)) THEN
   CALL CPU_TIME(time)
time=time-dummy
END IF


END SUBROUTINE DEFINE_TOPO


REAL (KIND=qPRec) FUNCTION gf(x,y,zeta,i,j)

REAL (KIND=qPRec), INTENT (IN) :: x,y,zeta
INTEGER, INTENT (IN) :: i,j

REAL (KIND=qPRec) :: h,hpx,hpxx,hpy,hpyy,z


h=hb(x,y)

hpx=(hb(x+deltax,y)-hb(x-deltax,y))/(2*deltax)

hpxx=(hb(x+deltax,y)-2*h+hb(x-deltax,y))/(deltax**2)

hpy=(hb(x,y+deltay)-hb(x,y-deltay))/(2*deltay)

hpyy=(hb(x,y+deltay)-2*h+hb(x,y-deltay))/(deltay**2)

z=zeta

SELECT CASE ((i-1)*3+j)
CASE(-3)! 00 special case, determinant of jacobian
gf=(dmax-h)/dmax
CASE(1) ! 11
gf=(1.-h/dmax)   
CASE(5) ! 22
gf=(1.-h/dmax)
CASE(9) ! 33
gf=(1./(1.-h/dmax)**2+(z-1)**2*(hpx**2+hpy**2)/(dmax-h)**2)*(1.-h/dmax)
CASE(2,4) ! 12 or 21
gf=zero   
CASE(3,7) ! 13 or 31
gf=1.*(z-1)/(dmax-hb(x,y))*(hb(x+deltax,y)-hb(x-deltax,y))/(2*deltax)*(1.-h/dmax)
CASE(6,8)  ! 23 or 32
gf=1.*(z-1)/(dmax-hb(x,y))*(hb(x,y+deltay)-hb(x,y-deltay))/(2*deltay)*(1.-h/dmax) 
END SELECT

END FUNCTION gf

ELEMENTAL REAL (KIND=qPrec) FUNCTION g(i,j,k,l,p)
INTEGER , PARAMETER :: alpha = 1
! 0<alpha <=1 controls the stretching of the mesh
! returns \sqrt(g)glp at location i,j,k
! NOTE THAT WHEN l==p IT IS RETURNED AT THE CORRESPONDING
! CELL EDGE
INTEGER, INTENT (IN) :: i,j,k,l,p

REAL (KIND=qPRec) :: xi,eta,zeta




SELECT CASE( (l-1)*3+p)
CASE (-3) !\sqrt(g)
xi=(i-.5)*deltax+xmin
eta=(j-.5)*deltay+ymin
zeta=(k-.5)*deltaz+zmin

g=alpha*(1.-h(i,j))*(abs(zeta))**(alpha-1)

CASE (1) !g11
xi=(i-1)*deltax+xmin
eta=(j-.5)*deltay+ymin
zeta=(k-.5)*deltaz+zmin

g=(1.-(h(i,j)+h(i-1,j))/(2))*(alpha*abs(zeta)**(alpha-1))
CASE (5) !g22
xi=(i-.5)*deltax+xmin
eta=(j-1)*deltay+ymin
zeta=(k-.5)*deltaz+zmin

g=(1.-(h(i,j)+h(i,j-1))/(2))*(alpha*abs(zeta)**(alpha-1))
CASE (9) ! g33
xi=(i-.5)*deltax+xmin
eta=(j-.5)*deltay+ymin
zeta=ABS((k-1)*deltaz+zmin)

g=(1.+(zeta**alpha-1)**2*(hx(i,j)**2+hy(i,j)**2))/(1.-h(i,j))&
     *(abs(zeta)**(1-alpha))/alpha



CASE (2,4) ! g12

xi=(i-.5)*deltax+xmin
eta=(j-.5)*deltay+ymin
zeta=(k-.5)*deltaz+zmin

g=zero


CASE (3,7) ! g13

xi=(i-.5)*deltax+xmin
eta=(j-.5)*deltay+ymin
zeta=ABS((k-.5)*deltaz+zmin)

g=(zeta**alpha-1)*hx(i,j)

CASE (6,8) ! g 23

xi=(i-.5)*deltax+xmin
eta=(j-.5)*deltay+ymin
zeta=ABS((k-.5)*deltaz+zmin)

g=(zeta**alpha-1)*hy(i,j)

END SELECT
END FUNCTION g

SUBROUTINE IND2CORD(Info,i,j,k,x,y,z,xi,eta,zeta)
! User definitions of data structures associated with each node
  ! Interface declarations
  TYPE (NodeInfo) :: Info        ! Data associated with this grid
  ! Internal declarations
  REAL (KIND=qPrec) hs
  integer i,j,k
  REAL (KIND=qPrec) x,y,z,xi,eta,zeta
  ! given indexes i,j,k, returns position in physical space and computational space
  
  xi=Info%xLower(1)+(i-0.5)*Info%dX(1)
  eta=Info%xLower(2)+(j-0.5)*Info%dX(2)
  zeta=Info%xLower(3)+(k-0.5)*Info%dX(3)
  
  x=xi*dmax
  y=eta*dmax
  z=((1-h(i,j))*zeta+h(i,j))*dmax
  
END SUBROUTINE IND2CORD

REAL (KIND=qPrec) FUNCTION hb(x,y)
REAL (KIND=qPrec) :: x,y

hb=.05*exp(-((x-15)/5)**4-((y-10)/2)**4)
!hb=.5*exp(-((x-10)/2)**4-((y-15)/5)**4)
!hb=.01*exp(-((x-15)/5)**4)
!hb=.01*exp(-((y-10)/2)**4)
!hb=.1*x/30
!hb=zero
END FUNCTION hb

LOGICAL FUNCTION g_mask(i,j)

INTEGER, INTENT (IN) :: i,j
SELECT CASE ((i-1)*3+j)
CASE(1) ! 11
g_mask=.true.   
CASE(5) ! 22
g_mask=.true.   
CASE(9) ! 33
g_mask=.true.   
CASE(2,4) ! 12 or 21
g_mask=.false.   
CASE(3,7) ! 13 or 31
g_mask=.true.   
CASE(6,8)  ! 23 or 32
g_mask=.true.   

END SELECT
END FUNCTION g_mask
END MODULE metric
